
function results = ocp(Params)
    import casadi.*

    % Create optimization problem
    opti = Opti();

    % Decision variables
    c  = opti.variable(Params.T,1);      % consumption
    ik = opti.variable(Params.T,1);      % capital investment
    is = opti.variable(Params.T,1);      % satellite investment
    k  = opti.variable(Params.T+1,1);    % capital stock
    s  = opti.variable(Params.T+1,1);    % satellite stock
    W  = opti.variable(Params.T+1,1);    % derelict satellites
    Z  = opti.variable(Params.T+1,1);    % rocket bodies
    F  = opti.variable(Params.T+1,1);    % fragments

    % Initial conditions
    opti.subject_to(k(1) == Params.k0);
    opti.subject_to(s(1) == Params.s0);
    opti.subject_to(W(1) == Params.W0);
    opti.subject_to(Z(1) == Params.Z0);
    opti.subject_to(F(1) == Params.F0);

    Utility = 0;

    for t = 1:Params.T
        
        % Output
        y = Params.a(t) * k(t)^Params.alpha1 * s(t)^Params.alpha2 * Params.N(t)^(1-Params.alpha1-Params.alpha2);

        % Resource constraint
        opti.subject_to(c(t) + ik(t) + is(t) == y);

        % Launches and collisions
        launches = Params.mu * is(t) / Params.q(t);
        x = Params.theta * s(t) * (W(t) + Z(t) + F(t));

        % Capital dynamics
        opti.subject_to(k(t+1) == (1-Params.deltak)*k(t) + ik(t));

        % Satellite dynamics (with collisions)
        opti.subject_to(s(t+1) == (1-Params.deltas)*s(t) + Params.q(t)*(1-Params.m(t))*is(t) - Params.theta*(W(t)+Z(t)+F(t))*s(t));

        % Derelict satellite dynamics
        W_next = (1-Params.deltaw - Params.epsilonw)*W(t) + Params.xi*Params.deltas*Params.mu *s(t)/Params.q(t) - Params.theta*(W(t)+Z(t)+F(t)+Params.mu*s(t)/Params.q(t))*W(t);
        opti.subject_to(W(t+1) == W_next);

        % Rocket bodies dynamics
        Z_next = (1-Params.deltaz-Params.epsilonz)*Z(t) + Params.psi*Params.mu*(1-Params.m(t))*is(t)/Params.eta - Params.theta*(W(t)+Z(t)+F(t)+Params.mu*s(t)/Params.q(t))*Z(t);
        opti.subject_to(Z(t+1) == Z_next);

        % Fragments dynamics
        F_next = (1-Params.deltaf)*F(t) ...
                 + (1+Params.nu1)*Params.omega*Params.mu*(1-Params.m(t))*is(t)/Params.eta ...
                 + (1+Params.nu1)*Params.phiw*Params.epsilonw*W(t) ...
                 + (1+Params.nu1)*Params.phiz*Params.epsilonz*Z(t) ...
                 + (1+Params.nu1)*Params.gammas*Params.theta*(W(t)+Z(t)+F(t))*Params.mu*s(t)/Params.q(t) ...
                 + (1+Params.nu1)*Params.gammaw*Params.theta*(W(t)+Z(t)+F(t))*W(t) ...
                 + (1+Params.nu1)*Params.gammaz*Params.theta*(W(t)+Z(t)+F(t))*Z(t);
        opti.subject_to(F(t+1) == F_next);

        % Bounds
        opti.subject_to(c(t) >= 1e-5);
        opti.subject_to(ik(t) >= 1e-5);
        opti.subject_to(is(t) >= 1e-5);
        opti.subject_to(k(t) >= 1e-1);
        opti.subject_to(s(t) >= 1e-1);
        opti.subject_to(W(t) >= 1e-1);
        opti.subject_to(Z(t) >= 1e-1);
        opti.subject_to(F(t) >= 1e-1);
        
        % Final State Bounds
        % opti.subject_to(k(end) >= 1e-10);
        % opti.subject_to(s(end) >= 1e-10);
        % opti.subject_to(W(end) >= 1e-10);
        % opti.subject_to(Z(end) >= 1e-10);
        % opti.subject_to(F(end) >= 1e-10);


        % Utility
        if Params.sigma == 1
            u = log(c(t)/Params.N(t));
        else
            u = (c(t)/Params.N(t))^(1-Params.sigma)/(1-Params.sigma);
        end
        Utility = Utility - Params.beta^(t-1) * Params.N(t) * u;
    end

    
    % Objective
    opti.minimize(Utility);

    % Solver options
    opts.ipopt.print_level = 5;
    opts.print_time = true;
    opti.solver('ipopt', opts);

    % Initial guesses
    opti.set_initial(c, Params.y0 * 0.8);
    opti.set_initial(ik, Params.k0 * Params.deltak);
    opti.set_initial(is, Params.s0 * Params.deltas);
    opti.set_initial(k, Params.k0);
    opti.set_initial(s, Params.s0);
    opti.set_initial(W, Params.W0);
    opti.set_initial(Z, Params.Z0);
    opti.set_initial(F, Params.F0);
       

    % Solve
    sol = opti.solve();

    % Output results
    results.c = sol.value(c);
    results.ik = sol.value(ik);
    results.is = sol.value(is);
    results.k = sol.value(k);
    results.s = sol.value(s);
    results.W = sol.value(W);
    results.Z = sol.value(Z);
    results.F = sol.value(F);
    results.utility = sol.value(Utility);
end
